Scaling of the thermal spectral function for quantum critical bosons in one dimension 
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We present an improved scheme for the precise evaluation of finite-temperature response func- 
tions of strongly correlated systems in the framework of the time-dependent density matrix renor- 
malization group. The maximum times that we can reach at finite temperatures T are typically 
increased by a factor of two, when compared against the earlier approaches. This novel scheme, 
complemented with linear prediction, allows us now to evaluate dynamic correlators for interacting 
bosons in one dimension. We demonstrate that the considered spectral function in the quantum 
critical regime with dynamic critical exponent z = 2 is captured by the universal scaling form 
S(k,ui) — 1/T ■ &s(k/VT,w/T) and calculate the scaling function precisely. 
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I. INTRODUCTION 

Response functions (B(t)A) quantify the effect of dis- 
tortions A of the system at time zero on the expec- 
tation values of observables B at time t. They con- 
tain important information on the governing quantum 
many-body physics [1] and are accessible in many differ- 
ent experimental setups. For example, recent advances 
in neutron-scattering techniques make very precise mea- 
surements possible [2] . It is of high importance to have 
numerical tools at hand that allow for an efficient and 
highly accurate computation of response functions for 
(strongly-correlated) condensed matter models in order 
to match theoretical models to actual materials and to 
gain an understanding of the underlying physical pro- 
cesses. Arguably, sufficiently precise tools are available 
for one-dimensional (ID) systems at zero temperature 
[3-6]. At finite temperatures, which, needless to say, are 
the relevant case for most experiments, our numerical 
abilities were however quite limited. As discussed and 
demonstrated in several works [7-9], finite-temperature 
response functions for strongly-correlated ID systems can 
be evaluated up to some maximum reachable time by us- 
ing the time-dependent density matrix renormalization 
group (tDMRG) [10-12]. A difficulty in the simulations 
of time-evolved states is the growth of entanglement with 
time [13-15]. In tDMRG calculations, this leads to a cor- 
responding severe increase of the computation cost and 
a strong limitation of the maximum reachable times, de- 
pending on the available computational resources. The 
effect is much more drastic for mixed states: At low tem- 
peratures, the corresponding computation cost basically 
increases by a power of two in comparison to the simula- 
tion of pure states. It is decisive to reach times that are 
sufficiently long to allow for the extraction of the desired 
physical information like spectral properties. 

In this paper, we present a novel tDMRG scheme for 
the calculation of dynamical correlators at finite temper- 
atures which typically doubles the maximum reachable 



times in comparison to the earlier schemes in the litera- 
ture. This allows for a very precise evaluation of thermal 
response functions which could not be addressed before. 
We also analyze and explain the computation cost for 
the different schemes. As a specific application we study 
bosons = Sij) with repulsive onsite and nearest- 

neighbor interactions as described by the extended Bose- 
Hubbard model 

H = -l^2{b\b i+1 +h.c)-^n i 

i i 

hih i+l . (1) 
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The new tDMRG scheme now allows us to examine the 
thermal spectral function [1] 
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with the partition function Zp — Tr e~^ H and the inverse 
temperature j3 = 1/T. The system becomes critical at 
the point fj, = — 1, T = as shown in Fig. 1. We will 
fix fj, = —1 in the following. In the z = 2 quantum- 
critical region of Fig. 1, for small quasi-momenta k and 
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Figure 1: Crossover phase diagram [16] around the quantum 
critical point at u = —1, T = 0. The T = state has zero 
density for fi < — 1, and non-zero density for fj, > — 1. In this 
paper, we compute the universal dynamic correlators in the 
z = 2 quantum critical region. Dynamic correlators in the 
dilute classical gas region were computed in Ref. [17]. 
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frequencies u>, the temperature T can be expected to be 
the dominating energy scale and we can expect 
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to be a function of and ^ 



i.e., S(k,u) « f(T) ■ 
<3> s (-^= ; As the integral of the spectral function over 
w (at fixed fc) is unity, we can conclude f(T) = A. 
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The existence and possible form of the scaling function 
$5 has been a long-standing open question. $s is ex- 
pected to be universal for certain (universality) classes 
of systems, as the long-ranged correlations at the criti- 
cal point should not depend on the specific form of the 
short-ranged interactions. In the language of renormal- 
ization group, this means that the renormalization flows 
of many other systems are governed by the same fixed 
point. This allows for a simplification of the calculations 
by going to the limit U — > oo, for which the system is 
restricted to have at most one boson per site. To assert 
the universality, we have introduced the nearest-neighbor 
density-density interaction in the model (1) and will show 
that the same scaling function applies for several different 
values of V. It should be stressed that the new tDMRG 
scheme, described in the following, is generally applica- 
ble for finite-temperature real-time and frequency-space 
simulations for strongly correlated ID quantum systems. 



II. METHOD 
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Figure 2: Maximum bond dimension max,M, occurring in the 
computation of the response function Xab f° r = A = V 
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chain length L = 128, U oo, V = 1, half filling (p = 1), 
and truncation weights ep = 10 -12 , e t = 10 -10 . The contour 
lines mark the values max;Mi = 2 6 , 2 7 , 2 s , 2 9 , 2 10 , 2 11 , i.e., 
the times that can be reached with those maximum bond 
dimensions. Schemes A and B reach much shorter times than 
the scheme of Eq. (9), optimized with respect to t' and /?'. 



purifications [34] . Due to the isomorphism oiH&H and 
B(H), the space of linear maps on the Hilbert space Ti, 
the MPS occurring in scheme A are in one-to-one relation 
with corresponding matrix product operators (MPO), i.e., 
operators of the form 
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In previous tDMRG applications [7-9], the thermal 
density matrix pp := e~@ H jZp was encoded by a cor- 
responding purification [18-21] which is a pure state 
\pp) € H (8 "Haux with an auxiliary system H aux = H = 
span{|cr)} such that Tr aux \pp){pp\ = pp. A matrix prod- 
uct state (MPS) [22-24] representation of the purification 
\pp) can be obtained by employing tDMRG for an imag- 
inary time-evolution starting from the infinite tempera- 
ture state (p oc 1, \p ) oc ^ W) ® k)) [7-9, 21]. Ex- 
pectation values then take the form Tr ppX = (pp\X\pp). 
The scheme of Ref. [7] to calculate thermal response func- 
tions 

Xj&W, t) ■■= Tr (pp B{t)A) - Tr (p p e i6t Be~ i6t A) (5) 

consists in obtaining first the MPS purification \pp) by 
tDMRG imaginary-time evolution, and to subsequently 
do a tDMRG real-time evolution to obtain \pp,t) :— 

e~ i6t \pp) as well as \App,t) := e~ i6t A\pp). With these, 
one computes the response function by evaluating the 
overlap Xab(Pi = {Pfit t\B\App, t). In the following, we 
refer to this procedure as scheme A. In this context, it is 
actually to some extent superfluous to think in terms of 



where the A i " i are Mj_i x Mj matrices with Mq = 
Ml = 1. The sizes Mi are also called bond dimensions. 
In a short-hand notation, scheme A can be denoted by 



Tr ([e-^' 2 e i6t ]B[e- ikt Ae-^/ 2 ] 



(7) 



The square brackets indicate which parts in this expres- 
sion are approximated as MPOs [Eq. (6)] (correspond- 
ing to the aforementioned MPS purifications \pp,t) and 
\App,t)) and are obtained via tDMRG. In this notation, 
the modified scheme B used in Ref. [9] reads 



Tr ( [e-' 3 "' 2 ] B [e- i6t Ae^ 6 ' 2 e i6t ] 



(8) 



In each step of the tDMRG, the evolved operators 
X = X(t) are approximated by an MPO with bond di- 
mensions Mi — Mi(j3,t) that are as small as possible 
for a given constraint on the desired precision of the ap- 
proximation [4]. This precision is in each step of the 
algorithm controlled by the so-called truncation weight 
e = (||X trunc — A"|| 2 /||A||2) 2 . Due to such truncations, 
the results of the different evaluation schemes like (7) and 
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(8) differ slightly from the exact Xab(P^) °f ^q. (5). It 
is essential to keep the errors e in the MPO truncations 
controlled and small at all times. 

Schemes (7) and (8) are typically far from optimal. 
One has a lot of freedom in designing a scheme that 
is as efficient as possible. With efficiency we mean 
that the occurring bond dimensions Mj, which deter- 
mine the computation cost, are as small as possible 
for given /3, t, and e. With a non-singular operator 
T, the most generic splitting involving two MPOs is 

i Tr ( [e i&t BT\ [f- 1 e- iAi Ae- pA ] ) . In principle, one 

would now like to optimize for T in order to minimize 
the computation cost and, thus, maximize the maximum 
reachable time. However, in the cases we studied, such 
optimizations turned out to be inefficient. The required 
cost scaled exponentially with the system size, even when 
restricting ourselves to optimize only with respect to uni- 
tary operators. Hence, we confine ourselves to study the 
computation cost of the less general class of schemes 

X [ e -iH(t-t')^ e -(0-0')H e iH(t-t')^ (Q) 

as a function of t' and j3' . Fig. 2 compares the evolu- 
tion of the occurring bond dimensions as a function of 
temperature and time for the three different schemes, 
Eqs. (7)-(9). In many cases, scheme B [Eq. (8)] has 
some advantage over scheme A. In Ref. [9] it was pointed 
out that the involved MPO [e- iAt Ae-^ 6 / 2 e i6t ] is time- 
independent for the simple case A = 1 , whereas the com- 
putation cost for \e~ lHt Ae~ / 2 ] , occurring in scheme 
A, can increase with time, even in this trivial case. More 
generally, the following argument applies for all operators 
A with finite spatial support: Typical condensed matter 
systems are quasi-local [25, 26], i.e., the spatial support of 
operators like e~ zHt Ae lHt , occurring in scheme B, grows 
only linearly with time. More precisely, outside a cer- 
tain space-time cone, the evolved operator acts almost 
like the identity and does hence not change the entangle- 
ment in that region. The preconditions are that all terms 
in the Hamiltonian are short-ranged and norm-bounded 
[25, 26]. Nevertheless, as Fig. 2 indicates, scheme A is 
often advantageous at very low temperatures, especially 
for non-critical systems. However, the scheme of Eq. (9) 
optimized for t' and j3' outperforms the earlier schemes 
substantially. A more detailed discussion and analysis 
of the different evolution schemes will be presented else- 
where [34]. As the example in Fig. 2 indicates, the maxi- 
mum reachable time is a slowly varying function of log f3 
and almost concave. Hence one can reach almost optimal 
results with scheme C 

^Tr([e i6tB e-? A Be- iAtB ] [ e - iAtA A e -* 6 e iAtA ]), 
Zp \ ) 

(10) 



1 1 1 


T=1/16, DMRG 
T=1/8, DMRG 

\ T=1/16, exact 

H T=1/8, exact 

\ T=1/4, exact 

"4 T=1/2, exact 

\ T=1 , exact - - ~ 


real part ' 

imaginary part — > 

i i i 


i i i 



-8 -6 -4 -2 2 4 6 8 



Figure 3: The exact time-dependent spectral function S(k = 
0, t) plotted as a function of t ■ T for fj, = — 1 and V = and 
compared with corresponding tDMRG data (scheme C). 



which does not require any optimization. After the 
imaginary time-evolution that yields [e~z H ], one runs 

two real-time tDMRG simulations to obtain MPOs 
^.fftBg-ffljjg-iHtfl] and \^~iHt A ^ e -^H e iHt A y With 

Eq. (10) one then obtains Xab(@^a + The accu ~ 
racy of the MPOs should be kept under control during 
the whole simulation, for example, as described above 
by bounding the truncation error. If this is done prop- 
erly, it is of minor importance what specific tA and 
tg = t — tA are chosen to evaluate x JbC^>*) ^ or a §i ven 
time t. For the typical case A = B\ the maximum reach- 
able times for tA and ts are equal, and the total max- 
imum reachable time with this scheme is then twice as 
large as the maximum time of scheme B. For all simu- 
lations in this article we used scheme C with a fourth 
order Suzuki- Trotter decomposition and a time step of 
size At = 1/8 in the tDMRG. The truncation weights 
were fixed to ep — 10~ 12 in the imaginary time-evolution 
and et = 10 -10 in the real-time evolution. 



III. EXACTLY SOLVABLE CASE 

For U — > oo and V — 0, the model (1) can be mapped 
to a system of free fermions by application of the Jordan- 

Wigner transformation Sj = rij=i(~l) CiCj ' &i an( i the 
resulting Hamiltonian — | 5^(4^+1 + h.c.) - 
can be diagonalized exactly. Due to Wick's theorem, all 
correlation functions are in this case determined by the 
single-particle Green's function [1]. The response func- 
tion (2), in particular, can be calculated by evaluating 
Pfaffian determinants of matrices that contain elements 
of the single-particle Green's function [27-29]. We also 
use this exactly solvable case to prove the high accuracy 
of the new generically applicable tDMRG scheme. 
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Figure 4: Spectral function (3) for fx = — 1 and V = 0, 4 and 
|, rescaled according to the scaling hypothesis (4). The top 
right plot shows S(k,oj) for ^7= = ~, all others for k — 0. 
The dashed lines around .5 = show, for the tDMRG data 
(scheme C), the contribution of the linear prediction, i.e., the 
second term in Eq. (11). The simulations were stopped when 
the computation cost per time step exceeded a certain value. 
Depending on the temperature, this occurred at maximum 
bond dimensions maxiAfi between 1000 and 4000. 

IV. LINEAR PREDICTION 

As the simulations yield the time-dependent spectral 
function only on a finite time interval [—t max ,t max ], de- 
fined by the maximum reachable time, a direct Fourier 
transformation to S(k,ui), Eq. (3), contains ringing arti- 
facts. To avoid them, one can use filter functions which, 
however, result in an artificial broadening. Instead, we 
use a linear prediction [30, 31] which basically fits a su- 
perposition of damped harmonic oscillations to the data. 
In the DMRG context, this technique was employed first 
in Ref. [5] for T = and in Ref. [7] for T > 0. We then 
Fourier transform the spectral function after extrapolat- 
ing it in this way for times \t\ > t ma . x ; 

S(k,w)= [ dte iujt S(k,t) + [ dte iut S LP (k,t), (11) 

where Si 1 p(k,t) is the result of the linear prediction. 
Needless to say, the quality of the linear prediction, and 
hence the precision of S(k, to), depend strongly on t maK . 

V. RESULTS FOR THE BOSONIC SPECTRAL 
FUNCTION 

In all simulations presented here, we used a system 
size of L = 128 and compared against larger systems 




Figure 5: Fitting the scaling function (4) for k — with the 
ansatz $s(0,w = f ) = ow/(l + buj 3 ) yields a = 0.649 and 
b = 0.826. 



to ensure that finite-size effects are negligible. Fig. 3 
displays the exact solution for V = and a very clear 
convergence of the properly rescaled S(k, t) curves. For 
a lattice size of L = 128, finite-size effects emerge for 
temperatures below T « 1/32 as small wiggles at larger- 
times (not shown). The displayed data basically show 
the behavior in the thermodynamic limit. Fig. 4 shows 
the spectral function in the frequency domain, confirming 
the scaling hypothesis (4) by the collapse of the properly 
rescaled curves for low temperatures, i.e., when plotting 
T ■ S(k 7 uj) as a function of cj/T and k/VT. The re- 
sults for V = \ and i confirm the universality of the 
scaling function $5. The way in which it is approached 
depends however on V. In the limit T — > the rescaled 
curves for different V should coincide. We would like to 
stress that in earlier calculations, based on the tDMRG 
schemes A and B, we were not able to reach times that 
would allow for a proper extraction of $s. Especially for 
low temperatures, the much smaller i max in those calcu- 
lations required a considerably larger contribution of 5lp 
to S(k,u>) [Eq. (11)] and, furthermore, resulted in rela- 
tively big errors of the linear prediction [7, 30, 31]. This 
caused considerable distortions of the curves in compar- 
ison to the quasi-exact results displayed in Fig. 4 which 
were obtained with the novel scheme C. Of course, also 
with the older schemes, one can always reach longer times 
by increasing the truncation weights and et, as this 
reduces the occurring bond dimensions Mi. But for val- 
ues greater than the ones chosen here, ep — 10 -12 and 
e t = 10~ 10 , the precision of the resulting data quickly 
deteriorates as we have checked by comparison against 
the exactly solvable case. A remarkably good fit of the 
scaling function = 0,Q = is given by ansatz 

ad)/(l + buj 3 ) with a = 0.649 and b = 0.826 as shown 
in Fig. 5. For nonzero k/VT, the scaling function differs 
from this simple ansatz at small Q but still decays as cj~ 2 
for large Li. 

In the continuum limit of the model (1), the large time 
and distance asymptotics of the real-space correlation 
function gfor) := T-^(b x+XQ (tpj with £ := 
and t := tT/2 can be evaluated analytically for U — » 00 
and V = in the framework of the Riemann-Hilbert 
problem [32, 33]. The corresponding formula given in 
section XVI. 9 of Ref. [33] splits into a factor Co (A) that 
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Figure 6: Comparison of T~ 1 ' 2 {b x + Xo (£)&£ ) against the ana- 
lytical formula for the asymptotics as derived by the Riemann- 
Hilbert problem formalism for V — [32, 33]. 

only depends on the ratio A := £/2r and terms that de- 
pend on A and r; 

X e ^f-oo dM|2r( M -A)| In ^ _|_ 0(t~ 1/2 )), 

with cj){n) := e y2^_\ ■ The formula for Co (A) involves 



a double integral that can not be evaluated easily, but 
the remaining terms are unproblematic. Fig. 6 shows, 
for space-time lines of constant A = £/2r, a compar- 
ison of the analytical formula for the asymptotics of 
r )/Co(A) against our numerical results for g. Both 
are consistent with each other. 



VI. CONCLUSION 

The presented novel tDMRG scheme (10) for the eval- 
uation of finite-temperature response functions outper- 
forms earlier approaches substantially and should hence 
be the method of choice for future applications. In this 
study, it allowed us to demonstrate that the thermal 
bosonic spectral function in the quantum critical regime 
with dynamic critical exponent z = 2 obeys a univer- 
sal scaling form and to obtain the corresponding scaling 
function. 
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